Breathing synchronization in interconnected networks 
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The harmony of an orchestra emerges from the individual effort of musicians towards mutual 
synchronization of their tempi. When the orchestra is split between two concert halls communicating 
via Internet, a time delay is imposed which might hinder synchronization. We describe this type 
of system as two interconnected networks of oscillators with a time delay and analyze its dynamics 
as a function of the couplings and communication lag. We discover a breathing synchronization 
regime, namely, for a wide range of parameters, two groups emerge in the orchestra within the same 
concert hall playing at different tempi. Each group has a mirror in the other hall, one group is in 
phase and the other in anti-phase with their mirrors. For strong couplings, a phase shift between 
halls might occur. The implications of our findings on several socio-technical and biological systems 
are discussed. 



PACS numbers: 05.45.Xt, 



>.75.-k, 



.60. aq 



Technology has furnished us with global connectivity 
changing the functioning of cooperative work, interna- 
tional business, and interpersonal relationships. Today, 
it is possible to distribute an orchestra over two concert 
halls in different continents. Fast Internet connections 
would provide the communication infrastructure to prop- 
erly combine the sounds [1]. As there is always a physi- 
cal limit speed to information transport, the communica- 
tion between sub-orchestras is subjected to a time delay 
due to the distance between halls. As we discuss here, 
this time delay migh pose a real challenge to the syn- 
chronizability among musicians, which is vital to a suc- 
cessful performance. When isolated, each sub-orchestra 
would naturally play in unison. However, in our model 
orchestra, when listening to the musicians in the other 
sub-orchestra, the musicians in the same hall can split 
into two groups, playing with different tempi, leading 
to breathing synchronization. Interestingly, the partner 
in the other hall will be always playing with the same 
tempo but either in phase or anti-phase, depending on 
their tempo. 

Understanding the consequences of a communication 
lag is also of major concern in other fields [2-4]. For 
example, the plasmodium Physarum polycephalum, an 
amoeba-like organism consisting of a network of tubu- 
lar structures for protoplasm flow, naturally shows peri- 
odic variations in its thickness, a necessary feature for its 
survival. A controlled setup has been prepared by Taka- 
matsu et al. where two regions of the same organism 
have been physically separated by a certain distance with 
the possibility of fine tuning the communication between 
them [5,6]. Depending on the coupling strength and time 
delay, the two regions have been shown to present phase 
and anti-phase synchronization of the oscillatory thick- 
ness. This is precisely what we find here for the orchestra 
in the regime of strong influence of other musicians in the 
same concert hall. As discussed in the final section, this 
biological system might be a prototype to experimentally 



evaluate the different regimes reported here. In what fol- 
lows, we focus on the example of the orchestra but our 
study might also have impact on several biological and 
tcchno-social systems as, for example, functional brain 
networks, living oscillators, or coupled power grids, as 
discussed at the end of this paper. 

When playing the same piece, musicians in an orches- 
tra try to synchronize their tempi, i.e., they try to play 
the notes at the same pace. Thus, for example, one vio- 
linist focuses simultaneously on the tempi of instruments 
in the same hall and on the corresponding violin in the 
other hall. All the other musicians act in the same way. 
The synchronizability of this setup can then be discussed 
in the framework of interdependent networks of oscilla- 
tors. Each concert hall is modeled by a network, where 
nodes are musicians and links represent the interaction 
between them. Additionally, each musician also estab- 
lishes a special inter-network coupling with one partner 
playing the same instrument in the other network (con- 
cert hall). Intra- and inter-network couplings have dif- 
ferent time scales: while the intra-network interactions 
can be considered instantaneous, the inter-network ones 
have a time delay that depends on the distance between 
concert halls. Recent geometrical studies of coupled net- 
works with intra- and inter-network links have revealed 
novel features never observed for isolated networks [7]. 
In particular, it has been shown that the overall robust- 
ness is reduced [8, 9] and the collapse of the system occurs 
through large cascades of failures [10, 11]. Dynamic prop- 
erties of coupled networks have also been studied [12-18], 
but the impact of a time delay on their synchronization 
is still an open issue, which we will address here. 

The Kuramoto model is the standard theoreti- 
cal framework for studying synchronizability of net- 
works [19-34]. For illustration purposes, here we stay 
with the example of an orchestra. A population 9 of n 
Kuramoto oscillators is considered to be mutually inter- 
acting, such as the musicians in one concert hall trying 



to keep the same tempo. We consider a random graph 
of average degree four. Each musician (oscillator) ie 6 
is described by a phase 0i(t), representing her/his cur- 
rent position, and a natural tempo Wj, corresponding to 
the pre-defined tempo of the music. Since all musicians 
are playing the same piece, we assume w, = ujq. The 
actual tempo of a musician is defined as the time deriva- 
tive of the phase, 9i (i) . To play the piece harmoniously, 
musicians try to synchronize their tempi. This interac- 
tion can be modeled in terms of the Kuramoto model 
as Oi — Ld + a Y^j=\ A® sin (6j — 0,), where the sum goes 
over all other musicians {i 7^ j), a is the coupling strength 
between them, and A® is the connectivity matrix such 
that Afj — 1 if musician i is influenced by j and zero 
otherwise. For simplicity, we assume that the musicians 
are all playing at the same unitary amplitude, so that 

the state of each musician can be described by a phasor 
e i9 i (t) i 

The collective performance in one concert hall, 
namely, the synchronization of its network, is char- 
acterized here by the complex order parameter 
re(i)e 4 *^' = ^ 5Z™ e l9j '(*' , where the sum goes over all 
musicians, ^(t) is the average phase, and the amplitude 
< |?*e(i)| < 1 measures the global coherence, i.e., how 
synchronized the musicians are. If rs(t) — 1 all mu- 
sicians play the same note at the same time, while low 
values of re imply that a significant fraction of musicians 
are out of phase. 

We introduce now a second population T, also of n os- 
cillators, representing the second concert hall. Within 
this network, musicians are coupled in the same way. 
Each j £ r is coupled with the corresponding partner 
i € 0, forming the inter-network couplings. In analogy to 
oscillators in 0, the motion of each oscillator is described 
by a phasor e nj ^\ of phase Jj(t). The inter-network 
coupling is subjected to a time delay r, corresponding to 
the time required for information to travel between con- 
cert halls [35]. Previous studies introduced time delay 
among oscillators of the same population [36, 37]. Here 
we consider the competition between an instantaneous 
intra-network and a delayed inter-network coupling. In 
a nutshell, the performance of each musician is described 
by, 
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FIG. 1. (color online) The interactions between a strongly 
delayed inter-network coupling and a weak intra-network cou- 
pling create two communities of different frequencies in steady 
state, a, Snapshot of populations at two different time steps 
(black dashed vertical lines in b)), for ojo = 1.0, r = 1.53, 
Tin — 0.01, and <jex = 0.5. The vertical position of each 
oscillator represents its phase, from — it to n, and the color 
represents the frequency at steady state. Superposition of 
these two communities leads to breathing synchronization. 
b, Time evolution of the order-parameter of populations Q 
(blue) and V (red) composed of n = 305 oscillators each with 
Lvo — 1.0, t — 1.53, and ctex = 0.5. Two scenarios of weak 
intra-network coupling are represented: am = 0.01 (continu- 
ous lines) and <7tn = 0.001 (dashed lines). 
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where the superscript t — r indicates the instant when 
the phases are calculated, and ctex and aw are the in- 
ter and intra-network couplings, respectively. In music, 
frequency is typically related with the note. However, 
here we assume that the properties of the note are not 
relevant. Instead, the usual natural frequency of the Ku- 
ramoto model (ui) corresponds to the tempo and, there- 



fore, frequency and tempo will be used as synonymous. 

We observe that for two interconnected networks of os- 
cillators with time delay, a weak intra-network coupling, 
and random initial distribution of phases, two frequency 
communities emerge within the same network, each syn- 
chronized with its mirror in a breathing mode, as shown 
in Fig. 1(a). In the figure, the color describes the fre- 
quency and the vertical position the phase. The fre- 
quency synchronization within groups occurs with phase 
locking. Interestingly, inter-network coupled pairs of 
nodes oscillate with the same frequency (same color) but 
might be either in phase or anti-phase (phase shift of tt). 
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FIG. 2. (color online) Scatter plot for the matrix of frequency pairs of intra-network neighboring oscillators, for 500 samples 
of n = 750, with ojo = 2.75, r = 1.53, o"ex = 1.5, and various <tin, namely, 0.4 (a), 0.8 (b), 1.2 (c), and 1.6 (d). The center 
of each symbol represents a pair of (wj,Wj), and colors and sizes correspond to their relative occurrence in the dataset. Blue 
empty circles correspond to results for am = 0, while the purple-yellow color scale is used for results with <tin > 0. 



Consequently, the presence of these two frequency groups 
affects the perception of the new global oscillatory state, 
which we call breathing synchronization. Figure 1(b) 
shows the time evolution of the order parameters re and 
rr for each population, quantifying this breathing be- 
havior. For each curve, the maximum corresponds to the 
instant at which both groups of frequencies are in phase, 
while the minimum to an anti-phase between groups in 
the same network. Additionally, since for one frequency 
there is a phase shift of it between inter-network pairs 
of nodes, the minimum in one network corresponds, nec- 
essarily, to the maximum in the other. In the context 
of the orchestra this implies that the communication lag 
in the interaction between concert halls generates a ten- 
dency to split musicians in each sub-orchestra into two 
groups playing at different tempi, hindering harmony. 
Cohesion within each community affects the amplitude 
of the breathing, as indicated by the order parameters 
for different values of <7in in Fig. 1(b). The weaker the 
intra-network coupling, the smaller is this amplitude. 

The observed breathing behavior is in deep contrast 
with what is expected for an isolated network (ctex = 0) . 
For isolated networks, the classical Kuramoto model 
is recovered, with frequency and phase synchronization 
emerging at a critical coupling am = a m . Above this 
threshold, a macroscopic fraction of oscillators is syn- 
chronized, all with the same frequency and phase. The 
value of (7j N increases with the variance of the natu- 
ral frequency distribution. Since here we consider the 
same natural frequency for all oscillators (w, = cl>o), 
<7j* n — » 0. The group of synchronized oscillators has 
frequency lu = loq [19] and the order parameter re(t) 
(or rr{t)) saturates in time at a non-zero steady-state 
value, which is a monotonically increasing function of 
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Interestingly, in the case of coupled net- 



works, and for sufficient inter-network couplings, none of 
the two frequencies is wq. 



To better understand the breathing synchronization, 
and in particular the emergence of frequency groups, let 
us consider the case of two coupled oscillators with time 
delay. The analytical solution obtained by Schuster and 
Wagner [35] for this problem indicates that, depending 
on the initial phase difference between oscillators (see 
Fig. SI in the Supplemental Material [38]), the pair can 
synchronize with different frequencies w, which are solu- 
tions of, 



lo = luq — ctex sin {lot) . 



(2) 



In spite of oscillating with the same frequency in the 
stationary state, the two oscillators might either be in 
phase, if cos(ojt) > 0, or anti-phase, otherwise. In the 
case of inter-connected networks, in the limit am = 0, 
the stationary state is expected to include all possible 
solutions of Eq. 2. Surprisingly, our results with a weak 
coupling reveal instead two frequency groups with phase 
locking. Nevertheless, the observed frequencies are con- 
sistent with the solution of Eq. 2. 

As we show next, when the internal coupling (am) is 
further increased, breathing synchronization is no longer 
stable and each network is synchronized, in one of two 
other synchronization regimes. To systematically study 
the dependence on am, we analyze the frequency cor- 
relation among intra-network neighbors i and j. Fig- 
ure 2 shows the scatter plots of the pair (oji,ujj) for dif- 
ferent values of intra-network coupling strengths. The 
limit ctin = is represented by the blue empty circles 
in all panels and the radius corresponds to the relative 
population of pairs when considering several samples. In 
this limit, the oscillators have all one of two possible fre- 
quencies, with four possible combinations of frequency 
pairs. From the relative size of the circles, we observe 
that the lowest frequency (w ~ 2.3 for ujq — 2.75 and 
t = 1.53) is the most populated one. As am is in- 
creased, the interaction among oscillators in the same 
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FIG. 3. (color online) Phase diagram for delayed coupled net- 
works. Parameter space of two coupling strengths <tex and 
ctin showing that the prevalence of one frequency over the 
other changes according to the coupling strengths. The color 
of each point represents the logarithm of the relative areas 
of the stable frequencies histogram, as represented in the top 
inset. The lower inset exhibits the phase boundaries for differ- 
ent time delays. The dominant mechanisms of each region are 
labeled accordingly: Breathing, Kuramoto [19], Competing, 
and Supernode states. 



network promotes their synchronization and, above a cer- 
tain intra-network coupling, each network does not split 
into two groups anymore. Instead of the breathing be- 
havior, the traditional synchronization is observed and 
the order parameter of each network saturates in time. 
As shown in Fig. 2(a), for ctin = 0.4 most nodes are 
synchronized with the lowest frequency and therefore a 
large percentage of the pairs are in the left-bottom cor- 
ner. Similarly to the Kuramoto model, in this competing 
state, oscillators synchronize at a unique stable frequency 
(cj rj 2.3), which is a solution of Eq. 2. As ctin is fur- 
ther increased (Fig. 2(b)-(d)), due to the strength of the 
intra-network coupling, each network tends to behave like 
a supernode and, depending on the initial conditions, one 
of two frequencies is obtained, which is again a solution 
of Eq. 2. Further analysis across samples (See Fig. S2 
in Supplemental Material [38]) also shows that the av- 
erage phase displacement between pairs of oscillators in 
different networks reaches A = 7r for intermediary val- 
ues of ctin, and decreases again once the supernodes are 
formed (See Fig. S2(a) in Supplemental Material [38]). 
For large ctin, the supernodes can be either in phase or 
anti-phase and, therefore, the average variance within a 
network has a value between zero and ir (See Fig. S2(b) 
in Supplemental Material [38]). 

To summarize the effect of several combinations of pa- 
rameters, we plot in Fig. 3 the phase diagram in the space 
of the two coupling strengths (o~m and <jex)- To iden- 



tify each regime, we compute the amount of oscillators 
with steady frequency below and above the mean value 
of possible frequencies (see Section Methods in the Sup- 
plemental Material [38]), A\ and A 2 , respectively, over 
different samples (see top inset of Fig. 3). The color 
map of the main plot of Fig. 3 shows the ratio of these 
quantities. While the blue area represents the domain 
of ctin and ctex combinations that leads to the smaller 
frequency, the shades in red represent the two regions 
where two frequencies can be achieved. Note that, the 
nature of the two synchronization regimes in red is dif- 
ferent. The one in the left (lower am) is characterized 
by the breathing behavior due to the presence of two fre- 
quency groups within each network. By constrast, in the 
supernode regime all nodes within a network are in phase 
locking, with the same frequency and, therefore, the or- 
der parameter is constant in time in the steady state. 
In the bottom inset, we show the phase boundaries for 
different time delays. From this, one can also see that 
the transition between regimes changes substantially for 
different time delays. In the Supplemental Material [38] 
we show that the transitions between regimes with one 
and two stable frequencies are abrupt. 

The presence of a time delay between two coupled net- 
works of oscillators poses a new challenge to the global 
control of the system. We have shown that the inter- 
play between coupling and delay leads to states of either 
a unique or two possible synchronized frequencies. We 
have found that, even with a weak intra-network cou- 
pling, oscillators within the same network split into two 
frequency groups. Each group has a mirror one in the 
other network oscillating at the same frequency. How- 
ever, depending on their frequency, a group can be either 
in phase or anti-phase with its mirror in the other net- 
work, resulting in breathing synchronization. Also, we 
show that an arbitrary increase of the intra-network cou- 
pling is not an option to achieve phase and frequency 
synchronization regardless of its initial conditions. In a 
certain region of the parameter space, the intra-network 
coupling promotes the formation of two supernodes (one 
per network), and two frequencies become stable. 

As previously mentioned, it is possible to prepare con- 
trolled experiments to evaluate the existence of these dif- 
ferent regimes in biological systems. Takamastu et al. [5] 
have shown that the distance and interaction strength 
between regions of a plasmodial slime mold can be fine 
tuned. This organism is a network of tubular structures 
with periodic variations in the thickness. In the exper- 
imental study, the focus was only on the regime where 
the intra-region interaction is much stronger than the 
inter-region one. Using the same methodology, it is pos- 
sible to control the intra-region interaction and study the 
different regimes described here. In particular, it would 
be of interest to observe oscillations with two different 
frequencies within the same region due to the communi- 
cation lag with the other region, resulting in breathing 



synchronization. 

Another example where synchronization in interdepen- 
dent networks certainly plays a relevant role is the human 
brain. Being a highly modular structure, its coherent op- 
eration must rely on the independence of different brain 
modules, which are functionally specialized, as well as 
on their efficient connection to ensure proper informa- 
tion transmission and processing. In a recent study [39], 
it was shown that the optimal integration of these mod- 
ules, which can be interpreted as complex networks made 
of intra-network couplings, is achieved through the addi- 
tion of long-range inter-network ties, therefore behaving 
globally as a small-world system. Moreover, their experi- 
mental observations are also consistent with the fact that 
these inter-network couplings should be spatially orga- 
nized in such a way as to maximize information transfer 
under wiring cost constraints [40]. To accomplish multi- 
sensory integration in this intricate architecture of neu- 
ronal firing-oscillators [41], however, information origi- 
nating from distinct sensory modalities (vision, audition, 
taction, etc.) must ultimately be processed in a synchro- 
nized way. This is typically the case when the processing 
of a visual signal influences the perception of an auditory 
stimulus and vice- versa [42, 43]. 
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METHODS 

Equation 1 has been numerically solved using a fourth order Runge-Kutta method with discrete time steps 5t = 
0.003. The stable frequencies were computed at t max — 100, using the difference between phases after one St step. 
The natural frequency has been chosen as ojq = 1.00 in Fig. 1 and luq = 2.75 for Figs. 2-5. Different values of wo do 
not affect qualitatively the results. The same values of 8t and t max were adopted for all simulations in this study. 

In Fig. 1, Panels a) and b) are based on one pair of undirected random networks of average degree four and 305 
nodes in each. Oscillators in this figure have been simulated for r = 1.53, o^x = 1-5- Panel a) is based on ctin = 0.01. 

In Fig. 2, Panels a)-d) contain the simultaneous representation of 500 pairs of random networks of 750 nodes. Color 
and size of each point represents the relative occurrence in all data. Oscillators in this figure have been simulated for 
r = 1.53 and ctex = 1-5. 

Fig. 3 is a schematic representation based on the average over 300 pairs of undirected random networks of average 
degree four and 500 nodes in each. The upper inset is a graphical representation of the histogram of all stable 
frequencies. The lower inset contains the same study for different delays, also averaged over 300 pairs of undirected 
random networks of average degree 4 and 500 nodes in each. A cutoff of u = 3.00, the midpoint of the stable 
frequencies for <tin = 0, was used to determine the areas A\ and A2- To avoid noise, we consider only frequencies 
with a relative occurrence of more than 10%. Oscillators in this figure have been simulated for r = 1.53 in the main 
Panel and several delays in the lower panel. 



STEADY STATE FREQUENCY ACHIEVED FOR VARYING INITIAL CONDITIONS 
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FIG. SI: For a combination of A , the initial phase displacement between inter-network oscillators, and w , we map 
the final frequency uj achieved in the case of a weak intra-network coupling. The color scale represents the variable 
lo — ujq . As for our simulations phases are initially randomly distributed, blue and red areas can be seen as the size 
of the basin of attraction of different frequency solutions. Data is an average of 500 sampels of n = 576 oscillators 
simulated with ctin = 0.01 and ctex = 0.50. Main and lateral panels show different time delays. 



DEPENDENCE OF THE PHASE DIFFERENCE AND FREQUENCY VARIANCE ON THE STRENGTH 

OF THE INTRA-NETWORK COUPLING. 
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FIG. S2: Each point in both plots represents an average over 500 samples of different network sizes with ujq — 2.75, 
r = 1.53, <7ex = 1-5. a, Phase difference A between pairs of inter-network connected oscillators, b, Average variance 
of the frequency of oscillators in the same network. Panels a) and b) are averages over 500 pairs of random networks 
of 250 nodes, 500 pairs of 500 nodes, and 500 pairs of 750 nodes. The standard deviations in all cases are smaller 
than the symbols. 



FREQUENCY COMMUNITIES IN TRANSITION REGIONS 

We analyze the transitions between synchronization regimes shown in the diagram of Fig. 3 in the main text. The 
transitions between regions of one and two stable frequencies are not smooth (See Fig. S3). We investigate three 
regions of the parameter space where transitions between states are expected. All are consistently abrupt, i.e., a 
small difference in the coupling strength triggers a bifurcation in the possible stable frequencies. Differences in the 
transition point can be explained in the context of finite size effects. A certain combination of coupling strengths and 
time delays leads to synchronization at a single stable frequency, but a small difference in any of the parameters can 
give rise to two stable frequencies, with the possibility of breathing synchronization, and a strong dependence of the 
final state on the initial conditions. 
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FIG. S3: For certain combinations of coupling strengths, one can observe abrupt transitions in which global 
synchronization is lost (Panels a and c) or recovered (Panel b). Points are averages over oscillators in each 
frequency community. Simulations were performed with 300 samples with t = 1.53 and wo = 2.75. In Panel a, we 
show results for am = 1.5, while Panels b and c are for ctex = 1-5. Panels a)-c) are averages over 300 pairs of 
random networks of 1000 nodes, 300 pairs of 1500 nodes, and 300 pairs of 1950 nodes. The standard deviation in all 
cases are smaller than the symbols. To reduce noise, we consider only frequencies with a relative occurrence of more 
than 10% to calculate averages. 



